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We investigate the aggregation kinetics of a simulated telechelic polymer gel. In the hybrid 
Molecular Dynamics (MD) / Monte Carlo (MC) algorithm, aggregates of associating end groups 
form and break according to MC rules, while the position of the polymers in space is dictated by 
MD. As a result, the aggregate sizes change every time step. In order to describe this aggregation 
process, we employ master equations. They define changes in the number of aggregates of a certain 
size in terms of reaction rates. These reaction rates indicate the likelihood that two aggregates 
combine to form a large one, or that a large aggregate splits into two smaller parts. The reaction 
rates are obtained from the simulations for a range of temperatures. 

Our results indicate that the rates are not only temperature dependent, but also a function of the 
sizes of the aggregates involved in the reaction. Using the measured rates, solutions to the master 
equations are shown to be stable and in agreement with the aggregate size distribution, as obtained 
directly from simulation data. Furthermore, we show how temperature induced variations in these 
rates give rise to the observed changes in the aggregate distribution that characterizes the sol-gel 
transition. 

PACS numbers: 61.41. +e, 82.35.-x, 83.80.Kn, 87.15. nr 



I. INTRODUCTION 

Applications of polymer gels range from traditional 
thickeners for paint |l| to novel precision gels for med- 
ical applications 0. In all cases, their ability to form 
reversible cross-links plays a crucial role since it allows 
the system to remodel the internal structure under the 
application of external stresses. Often these gels con- 
sist of triblock copolymers, containing a longhydrophilic 
spacer and short hydrophobic end groups [2|-|5||. I n an 
aqueous solution, these end groups aggregate. At high 
enough concentrations, aggregates form junction points 
in a gel network. At lower concentrations, the hydropho- 
bic groups still aggregate, but separate micelles form. 
The sol-gel transition characterizes the cross-over from 
fluid like to solid like behavior. Besides concentration 
changes both lowering the temperature or the application 
of external stress can also induce such a cross-over. A 
unique feature of the material is that the junction points 
are not permanent, instead the aggregate size - as well as 
its location - changes over time. 

Understanding the structure of the gel network as well 
as its dynamics due to formation and breakup of aggre- 
gating end blocks is crucial to the development of optimal 
designed materials. To this end, we have constructed 
a computer simulation of a toy system of associating 
telechelic polymers, in which end groups can form ag- 
gregates. Each polymer chain, 8 units in length, is mod- 
eled as a bead-spring chain molecule [f| . By means of a 
Monte Carlo step, we allow a chain end group to bind 
chemically with another end and hence form reversible 
junctions. Aggregates contain end groups that are con- 
nected to each other by these junctions. Their size varies 
and is set by many factors, such as the diffusion rate, the 



temperature of the system, and steric constraints due to 
the connection of each end group to a chain molecule. 
In previous work, we have calculated the aggregate size 
distribution for a range of temperatures @ and charac- 
terized the sol-gel transition. Moreover, we studied the 
network topology in detail using graph theoretical con- 
cepts Q. In this work, we report on the kinetics of the 
aggregation process and analyze the results in terms of 
the reaction rates in the master equations. 

Several other groups have studied the sol-gel transition 
by means of computer simulations. However, as far as we 
are aware, all of them have concentrated on the changes 
in structural properties instead of the kinetics underly- 
ing the reversible aggregation processes. For instance 
Bedrov et. al. Q performed simulations in which chain 
end groups were given a greater affinity to each other 
than to all other groups. They find similar aggregate 
size distributions as in our work, however since there are 
no specific interactions that bind end groups together, 
reaction rates cannot be determined in this model. The 
same is true for a recent study by Pandmanabhan et. al. 
in which they investigate semi flexible polymer gels [Io| . 

This paper is organized as follows. In section [TTJ wc 
give the background of the simulation model, which is 
followed by a discussion of the master equations and re- 
action rates in section Mil In section IIV1 we report on 
the measured rates and show that they are consistent 
with the aggregate distributions. Section [V] discusses 
how changes with temperature in the reaction rates are 
related to those in the overall dynamics of the systems 
that characterize the sol-gel transition. This is followed 
by a conclusion in section ["VTl 
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II. SIMULATION 

We employ a hybrid molecular dynamic (MD) / Monte 
Carlo (MC) simulation of telechelic polymers. The simu- 
lation uses a course-grained model introduced by Kremer 
and Grest [(|. Polymers are modeled as a string of beads. 
Beads interact with each other through a Lennard- Jones 
(LJ) potential 
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LJ 
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(1) 

for r < r c and zero otherwise. To speed up the calcula- 
tions the cut-off distance was taken to be r c = 2 1 / 6 and 
hence only the repulsive component of the potential is 
contributing. Beads connected by the chain structure are 
bonded through an anharmonic finitely extensible nonlin- 
ear elastic (FENE) potential 
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for rij < Rq and infinite otherwise. The parameters used 
were i?o = 1.5a and k = 30e<7~ 2 . It has been observed 
within the simulations Q that the average length of a 
FENE bond is 0.97a and that the chains cannot cross 
through each other. The positions of the beads are up- 
dated in MD simulations by integrating the equations of 
motion using a 5th-order Gear predictor-corrector algo- 
rithm with a time step of 0.005t. Length, energy, and 
time are expressed in terms of the parameters (a, e, r) of 
the LJ potential. 

The simulation system is comprised of 1000 telechelic 
polymer chains, each chain being eight beads in length. 
The simulation cell has dimensions 24 x 21 x 27 and hence 
the density is far below that at which the system turns 
glassy due to caging effects Q • Periodic boundary condi- 
tions are employed in the first two directions, while solid 
surfaces confine the simulation in the third direction. The 
beads at the chain ends arc functionalizcd linkers and can 
form reversible junctions with multiple other end groups. 
Formation and breaking of junctions takes place through 
MC moves [ll[. Junctions arc modeled by a FENE po- 
tential with the same parameters as described earlier, 
to which a constant negative constant U assoc is added. 
Every O.lr, an attempt is made to form and break junc- 
tions between all end groups. For the probability of a 
formation success a Metropolis algorithm is used which 
depends upon a Boltzmann factor exp(— AU /kT), where 
AU = U FENE + Uassoc is the energy difference between 
the old and the potentially new state. The magnitude of 
Uassoc affects the overall dynamics within the simulation 
and hence the temperature at which gelation occurs. As 
in previous work [7|, we choose U aS soc = — 22, which im- 
plies that the micelle transition temperature is T m = 0.5. 
Note that this algorithm employs a simplification of ne- 
glecting the energy barrier Eb, between the bonded and 
the unbonded states, that is present in an experimental 



system. It assumes a two-state system in which chain 
ends are either bound or unbound. As has been pointed 
out by Hoy et. al. 12j , the effect of Eb would have been 



to slow down the rate of formation and breakage of in- 
dividual junctions by a factor exp^—Es/kT). However 
the equilibrium constant, defined as the ratio of these 
rates, is independent of barrier height. Hence equilib- 
rium properties such as aggregate size distributions are 
not influenced by our simplification. It should be kept in 
mind that all rates at which aggregates form or break, 
reported later in this paper, are dependent on the cho- 
sen form of the interaction potential between the chain 
end groups (a simple FENE bond to which a constant 
negative association energy is added). They are also de- 
pendent on the fact that we attempt to break and from 
junctions every O.lr. In this study we take these set- 
tings as given and aim to relate the structural changes 
observed when cooling the system to measured changes in 
formation and breaking rates. Since most of our results 
depend only on the values of the ratio of rates and hence 
are independent of the barrier, we believe this study to 
be relevant to experimental systems as well. 

Finally we emphasize, that in our model there is no 
limitation to the number of junctions that an end group 
can form, hence detailed balance is observed within our 
work. Hoy et. al. [l2| adapted the model to perform a 
study involving only binary interactions, and therefore in 
that case strict detailed balance was lost. 

The simulation is analyzed through a range of temper- 
atures by thermally cooling from an initially high tem- 
perature. The temperature is controlled by coupling the 
simulation cell to a heat bath according to the fluctua- 
tion dissipation theorem as described by Kremer et. al. 
0. At each desired temperature, the system is allowed 
to equilibrate for at least 5000r prior to acquiring data. 
The spatial locations of each bead within the system, 
along with end group junction data is then gathered for 
a number of time steps necessary achieve statistically 
consistent results. Data at lower temperatures are ob- 
tained starting from a well-equilibrated configuration at 
T = 1.5. At this temperature, the system is first run with 
Uassoc = 0, so no junctions form. Given that each chain 
is only eight beads long, equilibrium is reached quickly. 
Junctions between chain ends are then introduced and 
the system is equilibrated until the number of junctions 
and the size distribution of the clusters of end groups 
fluctuates around its equilibrium value. The system is 
then slowly cooled at a rate of 2.500r per AT = —0.1, 
in order to reach the desired temperatures. Data ob- 
tained by cooling the system from different initial states 
at T = 1.5 are observed to be identical within statistical 
fluctuations. 

At the temperatures considered in this study the sys- 
tem is quite mobile. For instance the diffusion rate of 
aggregates equals 3 x 10~ 4 1/r at T = 0.45 and approx- 
imately 10~ 2 l/ratT=1.0. A single junction has an 
average lifetime of approximately 10 3 r at T = 0.45. The 
reported data have been obtained from runs that are at 
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least ten times the inverse diffusion rate of the system; 
hence the system samples many different configurations 
during this computer experiment. 

The aggregate size distribution is displayed in FigQ] p k 
is the probability of finding an end group in an aggregate 
of size k, where an aggregate contains only end groups 
that are connected to each other through the reversible 
junctions. It is calculated as: 
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(3) 



where N k is the average number of aggregates of size 
k and is obtained from the simulation data. N to t, the 
average total number of aggregates, depends on the tem- 
perature. 
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FIG. 1. (Color online) The time averaged aggregate prob- 
ability distribution of the polymer system for a range of 
temperatures. 

At high temperatures p k decreases monotonically with 
increasing cluster size. However, at T = 0.6 a shoul- 
der has formed, which develops in a pronounced peak 
at lower temperatures. The temperature at which the 
shoulder first forms (approximately T = 0.61) is called 
the onset temperature @, d, HH- At this temperature 
the properties of the system first start to deviate from 
those of an ordinary liquid. For instance, we have shown 
in previous work [7| that below this temperature the de- 
pendence of the relaxation time on temperature is not 
any longer Arrhcnius. When cooled even further the sys- 
tem undergoes a micelle transition. The total number of 
junctions between end beads <E> grows sharply and as a 
result a finite preferred aggregate size exists below the 
micelle transition. At this point, the distribution p k ex- 
hibits a distinct peak. The micelle transition tempera- 
ture is defined as the point at which d 2 $/dT 2 = fljj ]. 
The specific heat peaks as well. For the system at hand, 
we have calculated Q that the micelle transition tem- 
perature T m = 0.51. The micelle transition is an equi- 
librium transition in the systems structural properties. 



In addition at the gel transition the system becomes dy- 
namically arrested due to the long times that the end 
groups reside within specific aggregates. This nonequi- 
librium transition in the dynamics of the system occurs 
in our simulations at T = 0.4. 

In this paper, we will investigate by how far these tran- 
sitions in the properties of a gel forming system can be 
explained by the reaction rates at which aggregates of 
end groups combine and break up. 



III. KINETIC MODEL 

Reversible aggregation processes can be described as 
two competing effects: particle aggregation and cluster 
fragmentation. The traditional starting point fl6j for 
treating aggregation is a set of equations that describe 
how the sizes of the aggregated clusters change over time, 
called the master equations. These equations contain the 
rates at which aggregates form and break apart. In our 
simulations, an aggregate forms when two smaller ag- 
gregates closely approach each other and a new junction 
between chain end groups connects them. It breaks apart 
when all the junctions between two subunits of an aggre- 
gate disconnect. 

As an example of this process, we refer to a cartoon 
describing the breaking of an aggregate. At the center of 
Fig. six end groups (shown as red beads) are joined 
through multiple junctions. These end groups comprise 
an aggregate of size six. During the Monte Carlo step 
within the simulation, the two green junctions are bro- 
ken, resulting in two aggregates of size three (Fig. 03). 
The number of junctions broken, along with their distri- 
bution within the aggregate arc determined according to 
the rules of the Metropolis algorithm. Note that aggre- 
gation within the system is a reversible process, therefore 
the time reverse reaction will also occur. 
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FIG. 2. (Color online) A 2D cartoon illustrating the break- 
ing of an aggregate of size six (A) into two aggregates of size 
three (B). Note all simulations are preformed in 3D. 

If we denote an aggregate of size k by (fc), a change 
in the size of this aggregate can occur by the following 
reactions 



(t) + (k-£)^— f (k), 

q b ( k ,l) 

(k + £) " (k) + {£). 

qf (k+e,e) 



(4) 
(5) 
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In Eq. |4l qf(k,£) is defined as the rate that a given 
aggregate of size £ combines with one of size k — i to 
form an aggregate of size k within a time r. Similarly for 
the reverse reaction, £) is defined as the rate that an 
aggregate of size k breaks into an aggregate of size £ and 
one of size k — £. Of course an aggregate of size k is also 
formed when an aggregate of size k + £ breaks into one 
of size k and one of size £. The reverse reaction destroys 
an aggregate of size k. Eq. [5] describes the latter two 
events. 

The master equations describing the aggregation 
model can be written as 



dN k 



fc-i 



A— 1 



-£= ]T q f (k, £)N e N k _ e - Qb(k, £)N k 



+ ^q b (k + £,£)N k+e -J2q f (k + £,£)N k N e . (6) 

Here, N k is the number of aggregates of size k. The mas- 
ter equations represent a system of coupled differential 
equations. The model is finite since there are only 2000 
end groups in the simulated system. Due to steric effects 
the maximum observed aggregate size is even smaller. 
We never observed aggregates that contained more than 
approximately 60 end groups. 



IV. RESULTS 
A. Rates of Reactions 

The rates, at which the reactions according to Eq. 0] 
and Eq. [5] occur, are determined directly from the MD 
simulations. This is achieved by tracking the frequency 
of size specific reactions in both (fc) and (£). For in- 
stance, a formation occurrence Of(k,£) is defined as an 
event leading to an incremental change in the number of 
aggregates of size k in a period of time O.lr, caused by 
the combination of an aggregate of size £ and k — £. The 
time averaged reaction rates are then determined as the 
ratio of Of(k,£) to the number of aggregates which can 
produce (k). Hence, qf(k,£) is determined such that, 



q f (k,£) 



(7) 



The breaking rate qb(k,£) is determined using a similar 
procedure. 

Fig. [3] A, B, and C display color plots of the formation 
rates q f (k,£) at T = 0.45 (A), T = 0.55 (B), and T = 1.0 
(C). Rates shown in red are most reactive, where as dark 
blue is the least reactive. The breaking rates qb(k,£) 
are shown in (D), (E), and (F) for these temperatures 
respectively. Within the determination of the rates, it is 
assumed that (k) is always greater than (£) . It is observed 
that the reaction rates of the polymer system are not only 
temperature dependent, but are also a function of the size 
of aggregates involved in the reaction. 



The equilibrium constant, defined to be the ratio of 
the reaction rates Q(k 7 £) = qf(k,£)/qb(k,£), is related 
to the Gibbs free energy of the system and is shown at 
T = 0.45 (G), T = 0.55 (H), and T = 1.0 (I). Regions 
of large magnitude Q(k,£) mark (k) and {£) values, in 
which reactions arc more favorable to occur. 
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FIG. 3. (Color online) Average formation reaction rates 
q f (k,£) at a temperature of T = 0.45 (A), T = 0.55 (B), and 
T = 1.0 (C). The average breaking rates qb(k, £) are contained 
in (D), (E), and (F) respectively for each of the same temper- 
atures. The equilibrium constant Q(k,£) — qf(k,£)/qb(k,£) 
are displayed in (G), (H) and (I). Red indicates the most 
reactive (fc) and (£), while dark blue is the least reactive. 
Here, (£) is assumed never to be greater than (k). Average 
rates have been smoothed for display purposes. 



B. Consistency of the Master Equations with 
Aggregate Distribution 

The aggregate size distribution (Fig. [T]) should result 
from the master equations (Eqs. |6|) with the calculated 
reaction rates shown in Fig. [3] Two different methods 
are used to prove that this is indeed the case. The first 
consists of an unconstrained minimization of the master 
equations. At equilibrium dN k /dt = and the master 
equations become 



k-l 



k-l 



= J2lf( k ^) N e N k-e -^2q b (k,£)N k 
t=i e=i 

oo oo 

+ Y / 1b(k + £, £)N k+l - Qf(k + £, £)N k N e . (8) 



Eqs. [8] can be solved for iVt using Newton-based convex 
minimization techniques |17| . Using Eq. O p k is obtained 
from the solution for N k . The result at T = 0.55 is shown 
as stars in Fig. 2] Note that the solution is in excellent 
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agreement with the aggregate distribution obtained di- 
rectly from the data (circles). 




FIG. 4. (Color online) Solutions to the master equations 
of the kinetic model displaying agreement with the polymer 
aggregate distribution (circles) at T = 0.55. The dashed and 
dotted lines show the time evolution of the master equations. 
Starting from a constant (shown as dashed red), the distri- 
bution converges to the solid black line (R 2 = 0.999) after 
about 20t. Two intermediate steps are shown: the blue dot- 
ted line after lr and the green dot-dashed one after It. The 
triangles show that when only the {€) — 1 rates are employed 
the solution converges to the known aggregate distribution 
as well. However, the convergence is slower and occurs after 
500t. The convex optimization solution of Eqs. [8] is shown as 
stars. The inset is included so as to display the value of the 
solution at low k values. 

The second method consists of solving the master equa- 
tions in time. As an initial state, we take one in which 
each aggregate with k < 20 is equally likely and no ag- 
gregates with size over 20 exist. The initial distribution 
is shown as a dashed red line in Fig. 2) This distribu- 
tion is then updated by the master equations and the 
measured rates at T = 0.55 (the dotted blue line shows 
the distribution after time It, and the dot-dashed green 
line after time 2r) . The aggregate distribution converges 
in approximately 20r to that shown as a solid black line. 
As expected, this distribution agrees nicely with the mea- 
sured aggregate size distribution (shown as circles). Sim- 
ilar agreements were obtained for other initial distribu- 
tions and other temperatures. 



C. Simplification 

The condition of detailed balance [l3| can be used to 
show that only the rates of reactions, where one end 
group joins an aggregate or breaks away from it, com- 
pletely determine the aggregate distribution. This con- 
dition, if applicable, states that at a steady state the 
frequency of a specific reaction in the forward direction 
should be equivalent to that of the reverse reaction. In 



our case, that implies 

q f (k,£)N e N k _ e =q b {k,£)N k 
in which k 



for all (fc) and (£), 
follow that 

N, 



g/(M) 

q b (k,£) 



(9) 

> t. It would therefore 
N k - e N e . (10) 



This relationship at equilibrium has been verified within 
the simulations to hold (within statistical fluctuations) 
for multiple values of £, fc, and T. 

Eq. [TUJ has two implications in regards to the kinetic 
model. The first is that the first two terms and the last 
two terms in Eq. |8] add up to zero independently. The 
second is that the right hand side of Eq. [5] adds up to 
zero for each value of I independently. In other words, 
knowledge of q/(k, l)Ni/qb(k, 1) is sufficient to calculate 
the aggregate size distribution by propagating Eq. [101 for 
{£) = 1: N k = [q f {k,l)N 1 /q b (k,l)}N k _ 1 . The aggregate 
distribution p k can then be obtained from Eq. [3] Al- 
ternatively, starting from any initial condition, say the 
dashed red line in Fig. |4j and updating using only the 
(£) = 1 rates should yield the correct aggregate distribu- 
tion. This procedure was carried out and indeed yielded 
the correct distribution, see Fig. [4] (triangles). The con- 
vergence, however, is approximately a factor of 25 slower 
than the one attained when all rates are used. 

Although the (£) = 1 rates lead to the same correct 
distribution, the slow convergence here is due to the fact 
that they solely allow for a subset of the entire kinet- 
ics provided by the model using the complete range of 
rates. We make use of the fact that the (£) = 1 rates by 
themselves yield the observed aggregate distribution (see 
below) when we discuss how changes in these reaction 
rates give rise to the observed changes in the aggregate 
distribution that characterizes a micelle transition. 
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FIG. 5. (Color online) Magnitudes of the real component of 
the eigenvalues \i of the master equations' Jacobian matrix, 
evaluated at the known equilibrium state of the polymer at 
T = 0.55. The real eigenvalues are shown as red circles, while 
complex eigenvalues are shown in black squares. 
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D. Stability of Solutions 

To study the stability of the kinetic model under per- 
turbations near the actual solution, we analyze the eigen- 
values A,- of the Jacobian matrix [11 
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where yi represents the right hand side of dNi/dt ac- 
cording to Eq. [5] The size m of this matrix is deter- 
mined such that continuous, non-zero valued rates exist 
for 1 < k < m. Above that value of m, the rates are 
undetermined within the simulation time. This proce- 
dure results e.g. in m = 60 for T = 0.55. Since reaction 
rates and the aggregate size distribution are known, the 
Jacobian can be evaluated for a given temperature, and 
its eigenvalues calculated. 

The master equations consist of a multidimensional 
manifold. The resulting eigenvalues of Eq. [TT]are shown 
in Fig[5j They are either complex with negative real parts 
(black squares) or real and negative (red circles). This 
indicates that the solution is stable and that the singular 
point is a focus towards which the solutions are spiraling. 
Both the smallest and largest eigenvalues are real. The 
eigenvector corresponding to the largest absolute eigen- 
value peaks at N±. Hence, the fastest process consists 
of establishing the correct number of singles. The eigen- 
vector corresponding to the eigenvalue with the smallest 
absolute value is very similar to the equilibrium size dis- 
tribution. Apparently, the overall normalization of the 
solution is the slowest process [l9| . 



V. DISCUSSION 

A. Rates of Reactions 

Fig. OA displays q/(k,£) over a range of (k) values at 
T = 0.55. A symmetry inversion with respect to {£) is 
seen to exist within the formation reaction rates. For 
small (fc), the largest magnitudes of qf(k,£) are centered 
about (£) = (k)/2. This fact indicates, for (k) < 11, 
that two approximately equally sized aggregates partic- 
ipate in the formation reaction more frequently than an 
aggregate of small size £ and k — £. The behavior changes 
with increasing (k) and is displayed as an inversion of 
the symmetry within qf(k,£). For (k) > 11, where 
now the smallest magnitude of qf{k, £) is centered about 
{£) = (k)/2. Here, aggregates of small size £ associate 
with k — £ sized aggregates more frequently than two 
approximately equal sized aggregates. We label the (k) 
value, at which the symmetry inversion of <7/(fc, £) occurs 
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FIG. 6. (Color online) The average formation reaction rates 
qf(k,£) (A), breaking reaction rates qb(k, £) (B) as a function 
of (£) for several values of (k) at a temperature of T = 0.55. 
Errors within the data are on the order of the symbol size. 



by (fc)r- At T = 0.55, (k)r = 11 and increases with de- 
creasing temperature. The disassociation reaction rate 
<ft(fc, £), shown in Fig. [SJ3, does not exhibit such symme- 
try inversion. It is always easier to break off single beads 
than to split the system into equal parts. Note that the 
rates in figure [SK and |BJ3 are symmetric around (fc)/2 
since a split into aggregates of sizes {£) and (k — £) is 
identical to a split into sizes (k — £) and (£). 

Fig. [7] compares (&)t values for a range of tempera- 
tures with the aggregate size (k)s at which the aggre- 
gates are most spherical. The sphericity is calculated 
from the eigenvalues of the gyration tensor as described 
in [2(j. The smallest aggregates do not contain enough 
end groups to form a sphere. On the other hand, the 
largest aggregates have a long wormlike structure since 
steric effects caused by the polymer chains prevent the 
formation of a very large sphere. Thus, the most spheri- 
cal value appears at a middle aggregate size (fc)s- As can 
be seen both (fc) values in Fig. 0match, indicating that if 
an aggregate has a size larger than (fc)s, it is most likely 
to grow by addition of single beads or small aggregates 
at one of its ends. 



B. Dependence of Reaction Rates on Temperature 

Fig. [5] displays (Q) as a function of inverse tempera- 
ture. (Q) is calculated as a weighted average 



(Q) 



£ M Q(MMM) 



(12) 



where n(k,£) denotes the number of reactions involving 
aggregates of sizes k and £. Two separate cases are con- 
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FIG. 7. (Color online) The most spherically symmetric ag- 
gregate size (k)s (red circles), and the aggregate size where 
in qf(k,£) exhibits an inversion in symmetry (k)r (black 
squares) as a function of temperature. Lines are shown as 
guides for the eye. 



constant deviates from being described by an Arrhenius 
curve as the temperature decreases from T = 0.7. The 
low temperature data is fit by a stretched exponential 
(Q) = (Q)osxp(C/(T — To)) and is shown as a dashed 
line within the figure. The curve fitting results in pa- 
rameters T = 0.084, C = 2.78, and (Q) = 1.8 x 10~ 5 
for the {£) = 1 case, and T = 0.25, C = 1.49, and 
(Q)o = 1-3 x 10~ 4 for the all rates data. The extrap- 
olated temperature of T = 0.25 is in close agreement 
with our earlier published work Q, where the gelation 
temperature of To = 0.29 for a similar system was deter- 
mined from an analysis of relaxation times. The two fits 
cross at T = 0.65 for the {£) = 1 case and T = 0.6 when 
all (k) and (£) rates are used. This is close to the "onset" 
temperature, as was observed in previous work 

Again, we note that the (I) = 1 case describes a sub- 
set of the kinetics provided by the complete range of g's 
within the model. At T = 0.35 within the figure, the two 
cases converge. This is to be expected, since it has been 
observed that at this temperature, reactions are primar- 
ily composed of a single end group joining an aggregate 
or breaking away from it. 



sidered: (1) reactions in which only aggregates of size 
I = 1 are involved, and (2) all reactions. 
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FIG. 8. (Color online) A semi-log plot of (Q) as a function 
of 1/T, calculated for all (k) and (t) (red circles), and for 
(I) — 1 (black squares). The solid lines display an Arrhenius 
fit to the data within the high temperatures, while the dashed 
lines correspond to a stretched exponential fit within low tem- 
peratures. See the text for values of the fitting parameters. 



For high temperatures, Arrhenius behavior, as indi- 
cated by a linear trend on the semi-log scale, is ob- 
served. The rate data is fit by an Arrhenius curve 
(Q) = (<2)oexp(C/T) within the temperature range be- 
tween T = 1.2 and 0.7 and is shown as a solid line. In the 
case where (£) = 1 (shown as squares), the curve fit re- 
sults in parameters of C = 2.65 and (Q) = 4.49 x 10~ 5 , 
whereas when all (k) and {£) rates are used (shown as cir- 
cles), C = 3.26 and {Q) = 3.29 x 10~ 5 . The equilibrium 



C. Relation Between Reaction Rates and Size 
Distribution 

One of the signatures of a micelle transition is the oc- 
currence of a preferred finite aggregate size at low tem- 
perature [2l| characterized by the peak in the aggregate 
size distribution. We now investigate how this feature is 
related to the observed reaction rates. The only reactions 
considered are those in which one end group is added or 
removed from an aggregate. As shown above the rates for 
these reactions completely determine the aggregate size 
distribution. Fig. [5] shows the aggregate distributions 
together with these reaction rates for a range of temper- 
atures. Note that qb(k,l) increases faster with (fc) than 
qf(k, l)Ni. With decreasing temperature, qf(k, l)Ni be- 
comes less dependent on (k) and ultimately exhibits a 
decreasing trend with increasing (fc). This feature is nec- 
essary for the two rates to cross and provides a range of 
(fc), where in qf(k,l)N\ > qt,(k,l). For these (fc) val- 
ues, it is more favorable for a single end group to join 
the aggregate than to break from it. Moreover, at low 
temperatures, the minimum and maximum points in the 
aggregate size distribution are seen to correspond to the 
points where qf(k,l)Ni = 1). This is theoretically 
easily verified. An extremum in the aggregate distribu- 
tion implies that dp/dk — or pk ~ Pk~i- Using Eq. [31 
Eq. IHl can be rewritten as 



qf(k,£)N tot ptp k _ e = q b (k,£)p k , 



or 



Pk 

Pk-t 



q f (k,£) 
q b {k,£) 



q f (k,£) 
q b (k,£) 



N f 



(13) 
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FIG. 9. (Color online) The probability distribution pt and 
the crossing of the rates qf(k, l)N\ with qb(k, 1) at T = 1.0 
(A), T = 0.7 (B), T = 0.625 (C), T = 0.6 (D), T = 0.5 
(E), and T = 0.45 (F). At T = 0.625, the rates nearly cross. 
This temperature is in agreement with the onset of a shoulder 
within pk- For lower temperatures, the crossing of the rates 
is in agreement with the extrema of pk ■ The dashed lines act 
as guides for the eye and mark the (k) aggregate at which pu 
is at a maximum, in agreement with (k) when qf(k, l)JVi = 
q b (k, 1). 



for every {£). In particular, when i = 1 



Pk = g/(fc,l) 
~ q b [k,l) 



Pk-e 



N t ot Pi 



gKM) 
(ft (A, 1) 



(15) 



So pk = Pk-i when q/(k, 1)N\ — qb{k, 1). At high tem- 
perature q/(k, i)Ni < qb{k, 1) for all aggregates of size k 
and hence the aggregate size distribution has no extreme. 
The rates do not cross for T = 0.625 and higher, while 
they do for T = 0.6 or lower. The transition tempera- 
ture T ps 0.61 coincides with the temperature where the 
aggregate size distribution first forms a shoulder. This is 
expected since a necessary condition for the appearance 
of a peak in the aggregate size distribution is that the 
curves cross. Below T = 0.61, Arrhcnius dependence of 
the equilibrium constants and relaxation rates on temper- 
ature is lost. Moreover, the area between the qf(k, l)Ni 
and <?b(fc, 1) curves, shown in Fig[10j first increases and 
then decreases again when the temperature is lowered 
from T = 0.6. The area has an extremum at T = 0.5: 
the micelle transition temperature. 



Note that the above observation follows from the fact 
that both the forming and breaking rates are non mono- 
tonic as a function of (k). This fact deserves some fur- 
ther explanation. Most importantly, the two rates de- 
pend differently on (fc), one being concave up and the 
other increasing or concave down. We believe that these 
differences are due to cooperative effects. In general a 
single bead is attached to an aggregate with more than 
one junction (see Fig. [2]). In order for the bead to 
break away from the aggregate, multiple junctions have 
to break. The junctions have to break in a short time 
interval, since there is a significant chance that the first 
junction that broke would reform while the others are 
still there. On the other hand, when a bead joins an ag- 
gregate it is sufficient that one junction forms. qt,(k,l) 
initially decreases with (fc) since at small (fc) only one 
junction has to be broken and the need for two or more 
to be broken increases with (fc). In contrast, the larger 
the aggregate, the more places there are at which an end 
group can break off. Hence, at some point the trend re- 
verses and the rate starts to increase with (fc). On the 
other hand, 3/(fc, 1) initially increases with (k). This is 
due to the fact that a single bead can attach at more 
points to the aggregate of size (k — 1) as (k) increases. 
For this reason, one would expect qf(k,l) to increase 
monotonically. This is indeed the case at high tempera- 
ture. However at low temperature the rate starts to drop 
again at large (k) values. We believe that this is due to 
the much lower diffusion rate of large aggregates at low 
temperatures. 



x 10"' 




0.65 



FIG. 10. The area calculated as £ fe [qf(k, l)Ni - q b (k, 1)] 
over the range of (k) in which qf(k,l)Ni > qb(k, 1) as ob- 
served in Fig. 0(D) through (F) as function of temperature. 
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VI. CONCLUSION 

Usually in solution and other dense media one has to 
check the validity of transition rate theory, since if fric- 
tion is high, solvent dynamics cannot be neglected. This 
effect is very nicely demonstrated by Anna et.al. [22| . 
who showed experimentally how transition-rates change 
in different solvents, in good agreement with Kramers 
theory [23j . Kramers theory simulates solvent dynamics 
as frictional drag and influence of random forces. For the 
sake of simplicity we assumed here that both effects are 
small at the temperature range considered, and opted to 
use regular rate theory. 

We have investigated the kinetics of aggregation and 
fragmentation in a simulation of associating polymers. 
The process was described in terms of master equations 
and rates of reactions were obtained from the simulation. 
The master equations were shown to have a solution con- 
sistent with the aggregate size distribution as obtained 
directly from the simulations. Stability analysis showed 
that the aggregate distribution was indeed a stable so- 
lution to the master equations. We found that, since 
detailed balance holds within the simulation, knowledge 
of the {£) = 1 rates provides all the necessary information 
to reproduce the complete aggregate size distribution. 

It was observed that changes with temperature of the 
reaction rates can explain the changes in the aggregate 
distribution and dynamics near the sol-gel transition. 
The onset temperature of this transition was determined 



as the temperature at which there first exist a range of 
(k) values for which it is more likely that a single end 
group attaches to an aggregate than that one breaks off 
(qf(k,£)Ni > qb(k,£)). Near this temperature, (Q) be- 
gins to deviate from the Arrhenius, fluid-like behavior. 
In addition a shoulder develops within the distribution. 
Below the onset temperature, aggregates have a pre- 
ferred size. A symmetry inversion within the formation 
rate qj{k,£) occurs at this size and marks a fundamental 
change in the aggregate formation process. Aggregates 
smaller than the preferred size tend to be formed by join- 
ing two equal sized aggregates, whereas larger aggregates 
are formed by the addition of one end group. As tempera- 
ture continues to decrease, the likelihood that qf{k,£)N\ 
is larger than qbik, £) (measured as the area between the 
two curves of the rates versus (fc)) grows. The maximum 
at T = 0.5, is in agreement with the micelle transition 
temperature. 

In future work, we plan to investigate how the reaction 
rates are influenced by uniform and oscillatory shear. It 
is well known 24[ that the application of external stresses 
can influence the aggregate distribution; however its in- 
fluence on the kinetic processes in the system has been 
studied less. The master equations will be most likely a 
valuable tool in such studies as well. 
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